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ABSTRACT 

Evidence for an extraterrestrial flux of high-energy neutrinos has now been found in multiple searches 
with the IceCube detector. The first solid evidence was provided by a search for neutrino events with 
deposited energies > 30 TeV and interaction vertices inside the instrumented volume. Recent analyses 
suggest that the extraterrestrial flux extends to lower energies and is also visible with throughgoing, 
i/^-induced tracks from the Northern hemisphere. Here, we combine the results from six different 
IceCube searches for astrophysical neutrinos in a maximum-likelihood analysis. The combined event 
sample features high-statistics samples of shower-like and track-like events. The data are fit in up 
to three observables: energy, zenith angle and event topology. Assuming the astrophysical neutrino 
flux to be isotropic and to consist of equal flavors at Earth, the all-flavor spectrum with neutrino 
energies between 25 TeV and 2.8 PeV is well described by an unbroken power law with best-fit 
spectral index —2.50 ±0.09 and a flux at 100 TeV of ' 10“^® GeV“^s“^sr“^cm“^. Under the 

same assumptions, an unbroken power law with index —2 is disfavored with a significance of 3.8 a 
{p = 0.0066%) with respect to the best fit. This significance is reduced to 2.1 a {p = 1.7%) if instead 
we compare the best fit to a spectrum with index —2 that has an exponential cut-off at high energies. 
Allowing the electron neutrino flux to deviate from the other two flavors, we hnd a Vg fraction of 
0.18 ± 0.11 at Earth. The sole production of electron neutrinos, which would be characteristic of 
neutron-decay dominated sources, is rejected with a signihcance of 3.6 cr {p = 0.014%). 

Keywords: astroparticle physics — neutrinos — methods: data analysis 


1. INTRODUGTION 

With the completion of the IceCube detector in De¬ 
cember 2010, a decade-long journey towards building a 
kilometer-scale neutrino detector was finally concluded 
(for a review, see e.g. Spiering 2012). Using this in¬ 
strument, we aimed to detect the first extraterrestrial 
neutrinos, and thus to explore territories previously in¬ 
accessible to astronomy. These neutrinos have now been 
detected: In Aartsen et al. (2013, 2014e), the IceCube 
Collaboration presented the first evidence for a flux of 
high-energy extraterrestrial neutrinos, based on a search 
for events that start and deposit > 30 TeV inside the 
detection volume of IceCube. Indications for this flux 
had already been found in earlier analyses performed 
on data taken during the construction phase of IceCube 
(Aartsen et al. 2014c,d; Schonwald et al. 2013). More 


* Corresponding author (lars.mohrmann@desy.de). 


recently, the flux was measured down to deposited ener¬ 
gies of 10 TeV (Aartsen et al. 2015b), and evidence for it 
was observed using jz^-induced tracks from the Northern 
hemisphere (Aartsen et al. 2015a). 

Here, we present a maximum-likelihood analysis that 
is based on the event samples of all these analyses. The 
samples include track-like events, induced by charged- 
current interactions, as well as shower-like events, 
induced by charged-current Vg and Vr interactions and 
all-flavor neutral-current interactions.^ In this analysis, 
we derive improved constraints on the spectrum and the 
flavor composition of the astrophysical neutrino flux. 

In the remainder of this section, we summarize the 
phenomenology of astrophysical neutrinos and of back- 


^ In Ril7% of all charged-current Ut interactions, a muon is pro¬ 
duced in the tau decay, leading to a track-like instead of a shower¬ 
like signature. Furthermore, charged-current interactions may 
appear as shower-like if the muon leaves the detector unnoticed. 
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ground muons and neutrinos created in the atmosphere 
of the Earth, and we introduce the IceCube detector. 
The searches for astrophysical neutrinos that are used 
for the combined analysis presented here are discussed 
in section 2. The analysis method is explained in sec¬ 
tion 3, the results are given in section 4. We conclude 
with a discussion of the results in section 5. 

1 . 1 . Astrophysical Neutrinos 

Astrophysical neutrinos^ are created in interactions of 
high-energy cosmic rays with other massive particles or 
photons (Gaisser et al. 1995). The neutrinos, being elec¬ 
trically neutral and hence unaffected by cosmic magnetic 
fields, will travel in straight lines from their point of ori¬ 
gin to Earth. If the interactions happen close to the ac¬ 
celeration sites of the cosmic rays, they will thus reveal 
those. Unlike gamma rays, which are created in the same 
processes, the neutrinos are also unlikely to be absorbed 
during their journey (see e.g. Learned & Mannheim 
2000). Because of these properties, neutrinos are ideal 
messengers to study the sources of high-energy cosmic 
rays. Astrophysical neutrinos carry information abont 
these sources in their energy spectrum and flavor compo¬ 
sition, even if their individual positions on the sky cannot 
be resolved yet (Hooper et al. 2003; Choubey & Rode- 
johann 2009; Lipari et al. 2007; Laha et al. 2013). 
Candidate sources include active galactic nuclei (e.g. 
Stecker et al. 1991; Miicke et al. 2003), gamma-ray bursts 
(e.g. Waxman & Bahcall 1997; Guetta et al. 2004), star- 
burst galaxies (e.g. Loeb & Waxman 2006), and galaxy 
clusters (e.g. Murase et al. 2008), as well as galactic 
objects like supernova remnants or pulsar wind nebu¬ 
lae (e.g. Bednarek et al. 2005; Kistler & Beacom 2006; 
Kappes et al. 2007). 

To first order, the energy spectrum of astrophysical 
neutrinos follows that of the cosmic rays at their acceler¬ 
ation sites. If Fermi shock acceleration is the responsible 
mechanism, a power law spectrum E~'^ with 7 ~ 2 is 
expected (Gaisser 1990), although the details depend on 
the characteristics of the specific sources (see e.g. Becker 
2008). Furthermore, the majority of the astrophysical 
neutrinos are expected to arise from the decay of pions 
created in cosmic-ray interactions, i.e. tt —>■ fi fol¬ 
lowed by /r —)■ e + Vg + Vfj^. The flavor composition result¬ 
ing from this decay chain is r'e : = 1 : 2 : 0. Tak¬ 

ing into account long-baseline neutrino oscillations, the 
flavor composition at Earth is different, approximately 
1 : 1 : 1 for this scenario (Learned & Pakvasa 1995; 
Athar et al. 2006). This first order model of the energy 
spectrum and flavor composition has often been used as 
a benchmark scenario in the past. 

Second order corrections to this benchmark model arise 
e.g. from muon energy losses (Kashti & Waxman 2005) 
and muon acceleration (Klein et al. 2013), these effects 
can alter both the energy spectrum and flavor compo¬ 
sition of the astrophysical nentrino flux. While the en¬ 
ergy spectrum is difficult to constrain by general argu¬ 
ments, there are two limiting scenarios for the flavor 
composition at the sources^: muon-damped sources, in 


^ Here and in the rest of this article, we imply also anti-neutrinos 
when we speak of neutrinos. 

^ Neglecting the production of tau neutrinos at the sources, 
which is a common assumption (Choubey &; Rodejohann 2009). 


which the high-energy neutrino flux is suppressed due to 
energy loss processes of the muons, and neutron-beam 
sources, in which the neutrinos are created from neu¬ 
tron rather than pion decays. The flavor compositions 
at the source for these (idealized) scenarios are 0 : 1:0 
and 1:0:0, respectively (Lipari et al. 2007). Us¬ 
ing the neutrino oscillation parameters from Gonzalez- 
Garcia et al. (2014, inverted hierarchy), the expected 
flavor transitions for the three source classes discussed 
here are 1 : 2 : 0 —0.93 : 1.05 : 1.02 (pion-decay), 
0:1:0—)' 0.19 : 0.43 : 0.38 (muon-damped), and 
1:0:0 —) 0.55 : 0.19 : 0.26 (neutron-beam). 

The implications of deviations from the benchmark 
scenario have been discussed by several authors. As¬ 
suming extragalactic sources. Winter (2013) argues that 
a spectral index 7 > 2.3 is more easily explained in 
photohadronic scenarios, i.e. by py-interactions as the 
origin of the neutrinos. Similarly, Murase et al. (2013) 
point out that if a hadronuclear origin (pp-interactions) 
in extragalactic sources is assumed, measurements of the 
diffuse extragalactic gamma-ray background (see Ack- 
ermann et al. 2015) imply 7 < 2.1 — 2.2 under certain 
conditions. However, hadronuclear models that take into 
account the diffuse gamma-ray background as well as re¬ 
cent IceGube data (Aartsen et al. 2015b), which point 
towards a softer spectrum, have also been proposed, 
e.g. by Senno et al. (2015). Other authors invoke spec¬ 
tral arguments to propose that the astrophysical neu¬ 
trino flux is produced by sources within the Milky Way 
(Neronov & Semikoz 2014; Gaggero et al. 2015). Finally, 
the implications of flavor ratios different from the bench¬ 
mark scenario have been discussed e.g. by Beacom et al. 
(2003); Lipari et al. (2007); Vissani et al. (2013). Bus¬ 
tamante et al. (2015) give an overview over the flavor 
compositions resulting from various standard and non¬ 
standard neutrino production and propagation scenarios. 

1 . 2 . Atmospheric Backgrounds 

All relevant backgrounds to searches for high-energy 
astrophysical neutrinos are created in cosmic ray-induced 
air showers in the atmosphere of the Earth. Of all the 
particles created in air showers, only muons and neutri¬ 
nos can reach the IceGube detector. Atmospheric muons 
constitute, by far, the most abundant background, trig¬ 
gering the detector at a rate of several kHz. Gharacter- 
istically, they reach the detector from above the horizon 
and are first detected on the boundary of the instru¬ 
mented volume. Atmospheric neutrinos are created at a 
similar rate, but are detected much less frequently due to 
their small interaction probabilities. They reach the de¬ 
tector from all directions, in particular also from below 
the horizon. Atmospheric neutrinos that have passed 
through the Earth are free of muon background, but 
also hard to distinguish from astrophysical neutrinos. In 
contrast, atmospheric neutrinos arriving from above the 
horizon are often accompanied by atmospheric muons 
(Gaisser et al. 2014), which is never the case for astro- 
physical neutrinos. 

At lower energies, the atmospheric neutrino flux is 
dominated by so-called conventional atmospheric neu¬ 
trinos from the decays of kaons and charged pions. Since 
these particles are likely to interact with air molecules be¬ 
fore they decay, the resulting neutrino flux differs from 
the original cosmic ray flux in its energy and zenith an- 
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gle dependence: the energy spectrum is steeper (approxi¬ 
mately and the flux is enhanced towards the hori¬ 

zon (Gaisser 1990). An additional contribution to the 
atmospheric neutrino flux is expected from the decays 
of heavy, short-lived hadrons containing a charm or bot¬ 
tom quark. At the energies relevant here, these almost 
always decay before having the chance to interact, giving 
rise to a flux of prompt atmospheric neutrinos. This flux 
is predicted to follow that of the cosmic rays more closely, 
with an energy spectrum of approximately and an 

isotropic zenith angle distribution. It is, however, yet 
to be conclusively observed and predictions for its mag¬ 
nitude vary strongly (Bugaev et al. 1989; Martin et al. 
2003; Enberg et al. 2008; Bhattacharya et al. 2015). 

1.3. The IceCube Detector 

The IceCube detector (Achterberg et al. 2006) consists 
of 5160 optical modules, installed in the ice underneath 
the geographic South Pole at depths between 1450 and 
2450 meters. The modules are deployed on 86 strings, 
where the vertical spacing between two modules on a 
string is approximately 17 meters and the horizontal 
spacing between two strings is approximately 125 me¬ 
ters. Eight of the strings are part of a more densely 
instrumented region in the center of the detector, Deep- 
Core, which serves as a low energy extension to IceCube 
(Abbasi et al. 2012). Each optical module is a spherical 
glass housing that contains a 10"-photomultiplier (Ab¬ 
basi et al. 2010), together with digitization electronics 
(Abbasi et al. 2009). The IceCube Neutrino Observa¬ 
tory also includes a cosmic-ray detector called IceTop, 
which consists of 81 stations that are installed on the 
surface above the IceCube detector (Abbasi et al. 2013). 

IceCube detects neutrino interactions by measuring the 
Cherenkov radiation that is induced by the secondary 
particles created in the interaction. The signature of a 
neutrino interaction in the detector depends on the fla¬ 
vor of the incoming neutrino as well as on the interaction 
type. So-called tracks arise from charged-current in¬ 
teractions, in which a muon is produced together with 
a hadronic particle shower at the interaction vertex. At 
the relevant energies, the muon has a range of several 
kilometers, leading to an elongated, track-like signature. 
Due to the long lever arm, the directional reconstruction 
of these events is very precise, with a typical median reso¬ 
lution of better than 1° (Aartsen et al. 2014a). Since the 
muon will typically leave the detector, and often also en¬ 
ter it from outside, only a lower limit can be determined 
for the neutrino energy. On the other hand, charged- 
current Vg and i/r interactions as well as neutral cur¬ 
rent interactions of all flavors give rise to showers.'^ The 
electromagnetic and hadronic particle showers produced 
in these interactions have dimensions that are typically 
much smaller than the detector spacing, so they appear 
as point sources of light. As a result, the directional 
reconstruction of these events is worse than for tracks, 
typically 15° (Aartsen et al. 2014b). For both event 
topologies, the energy deposited in the detector, which 


^ At very high energies (> 1 PeV), the tau created in a charged- 
current Vt interaction can travel sufficiently far (~ 50 m/PeV) 
before decaying such that two separated showers may be observed; 
one at the interaction vertex and one at the decay point of the tau. 
However, this “double-bang” signature has not been observed yet. 


is a lower bound for the neutrino energy, can be recon¬ 
structed with a resolution of about 15% (Aartsen et al. 
2014b). For a track event, usually the differential en¬ 
ergy loss or the total muon energy are estimated. These 
quantities are correlated with the energy of the neutrino 
that produces the muon. However, their absolute values 
should be interpreted with care (Aartsen et al. 2014b). 
Because the absolute values are unimportant here, the 
quantities are quoted in arbitrary units in this article. 

The construction of the IceCube detector was finished 
in December 2010. However, data had already been 
taken since 2005, after the installation of the first Ice- 
Cube string. Data used in this analysis were taken with 
configurations consisting of 40 strings (2008/2009), 59 
strings (2009/2010), 79 strings (2010/2011) and the full 
86-string configuration (since 2011). 

2. SEARCHES FOR HIGH-ENERGY NEUTRINOS WITH 
ICECUBE 

Searches for high-energy astrophysical neutrinos have 
been performed ever since the beginning of IceCube data- 
taking, as before with data collected by IceCube’s pre¬ 
decessor AMANDA (see e.g. Achterberg et al. 2007; Ab¬ 
basi et al. 2011). Here, we consider only more recent 
searches, and only those targeted at the TeV-PeV en¬ 
ergy range. We also restrict ourselves to searches for a 
diffuse flux of neutrinos, i.e. we do not include searches 
that solely aim to identify individual sources, since these 
typically contain a larger fraction of atmospheric muon 
background events. This leaves a total of six searches, 
listed in Table 1. 

The searches mainly differ in their strategy to sup¬ 
press the large background of muons and neutrinos that 
are created in cosmic ray-induced air showers in the at¬ 
mosphere of the Earth. The important characteristics 
of such a strategy are the event topologies that are se¬ 
lected (listed in the table under “Topology”), the use 
(or lack) of containment criteria (“Containment”), and 
the resulting neutrino energy and zenith angle ranges in 
which the search is sensitive to a flux of astrophysical 
neutrinos (“Energy range”, “Zenith range”). The ta¬ 
ble also lists the corresponding data taking period and 
the observables that are used in the maximum-likelihood 
analysis. The individual searches and their background 
suppression strategies are briefly explained in the follow¬ 
ing sections, where we have classified them by topology 
into searches for track-like events (Tl, T2), searches for 
shower-like events (SI, S2) and hybrid event searches, 
which select both kinds of events (HI, H2). For more 
detailed information, we refer to the references listed in 
Table 1. A table listing the fraction of neutrino interac¬ 
tion types that contribute to the individual samples can 
be found in Appendix A. 

2.1. Searches for Track-Like Events 

Searches of this kind select upward-going and hori¬ 
zontal zz^-induced tracks. The good angular resolution 
of these events provides excellent discrimination against 
downward-going muon events, which can be reduced to 
a negligible level. Because events with interaction ver¬ 
tices outside the instrumented volume can be selected, 
this strategy allows the selection of large neutrino sam¬ 
ples of high purity. It is however restricted by design to 
directions from which only neutrinos can reach the detec- 
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Table 1 

Searches Combined in the Maximum-likelihood Analysis 


ID 

Topology 

Containment 

Energy range^ 
(TeV) 

Zenith range 
(deg) 

Data taking 
period 

Observables 

Reference 

Tl 

tracks 

no 

> 100 

90 - 180 

2009 - 2010 

energy, zenith 

1 

T2 

tracks 

no 

> 100 

85 - 180 

2010 - 2012 

energy, zenith 

2 

SI 

showers 

yes 

> 100 

0 - 180 

2008 - 2009 

energy 

o'" 

S2 

showers 

yes 

> 30 

0 - 180 

2009 - 2010 

energy 

4“ 

HI 

showers, tracks 

yes 

> 50 

0 - 180 

2010 - 2013 

energy, zenith 

5,6 

H2 

showers, tracks 

yes 

> 20 

0 - 180 

2010 - 2012 

energy, zenith, topology 

7 


References. — (1) Aartsen et al. (2014c); (2) Aartsen et al. (2015a); (3) Aartsen et al. (2014d); (4) Schonwald et al. (2013); 

(5) Aartsen et al. (2013); (6) Aartsen et al. (2014e); (7) Aartsen et al. (2015b) 

^ Refers to the approximate neutrino energy range in which the analysis is sensitive to astrophysical neutrinos, not the range of deposited 
energies selected. 

Samples “la” and “Ib” from the reference are used here. 

“ “Analysis A” in the reference. Also see Appendix B for additional changes to the event selection made within this work. 


tor, i.e. the northern hemisphere and directions close to 
the horizon. This implies that atmospheric neutrinos are 
an irreducible background to searches of this type and 
can only be distinguished from astrophysical neutrinos 
on a statistical basis. Because atmospheric neutrinos are 
predominantly muon neutrinos at all relevant energies, a 
flux of astrophysical neutrinos will only be visible at very 
high energies > 100 TeV. Neutrinos of these energies have 
a non-vanishing probability to be absorbed during their 
passage through the Earth (Gandhi et al. 1996), searches 
for track-like events are therefore most sensitive in the 
direction of the horizon, where the amount of traversed 
matter is less compared to more vertical directions. 

Here we use the neutrino event samples from two dif¬ 
ferent event selections following this approach, performed 
on data taken in different periods. The first one (“Tl”) 
uses data taken in 2009/2010 with the 59-string configu¬ 
ration of IceCube while the second one (“T2”) uses two 
years of data, taken between 2010 and 2012 with the 79- 
string and the 86 -string configuration of IceCube. For 
both searches, energy and zenith angle information are 
available for the neutrino sample. Note that we do not 
use the full event sample provided by the analyses, but 
restrict ourselves to the high-energy tail of the spectrum. 
This degrades the ability to constrain the atmospheric 
neutrino spectrum (and thus, indirectly, the astrophys¬ 
ical spectrum) somewhat, but facilitates the consistent 
treatment of systematic uncertainties (which are diffi¬ 
cult to account for in the high-statistics threshold region 
of the track-like searches) in the combined analysis pre¬ 
sented here. 

2.2. Searches for Shower-Like Events 

These searches select shower-like events and are thus 
sensitive to all neutrino flavors. Lacking the good an¬ 
gular resolution of tracks, searches of this type achieve 
background suppression by selecting events with a spher¬ 
ical topology, by requiring the interaction vertex to be 
contained within the instrumented volume of IceCube, 
and by vetoing events in which light is first detected 
on the outermost layer of optical modules. This means 
that the sensitive analysis volume is smaller than the in¬ 
strumented volume, but allows the selection of neutrinos 
coming from any direction. However, reliable rejection of 
muon background events is only possible if there is a suffi¬ 
ciently high probability of detecting the incoming muon. 
This typically leads to a contamination of residual muon 


background events at low energies, where muons may 
produce too little light to be rejected efficiently enough. 

Event samples selected by two searches for shower-like 
events are included here, one performed on data taken in 
2008/2009 with the 40-string configuration of IceCube 
(“SI”) and one performed on data taken in 2009/2010 
with the 59-string conhguration of IceCube (“S2”). Note 
that the event sample of search S 2 has been extended 
to lower energies compared to what was presented in 
Schonwald et al. (2013), see Appendix B for further de¬ 
tails. Neither analysis originally accounted for the at¬ 
mospheric self-veto effect, i.e. the possibility to reject 
atmospheric neutrinos due to atmospheric muons that 
accompany them (Schonert et al. 2009). For the anal¬ 
ysis presented here, the results were corrected for this 
effect based on a calculation by Gaisser et al. (2014); 
checks were made to ensure that systematic errors in the 
description of the effect for these event samples do not 
affect the results obtained here. The search SI provides 
two event samples with slightly different selection cri¬ 
teria, a low-energy sample (with deposited energies be¬ 
tween 2 TeV and 100 TeV, “Sla”) and a high-energy sam¬ 
ple (with deposited energies > 100 TeV, “Sib”). Both 
shower searches provide reconstructed energies for their 
respective event samples. 

2.3. Hybrid Event Searches 

The two hybrid searches considered here do not make 
any requirements on the topology of the selected events, 
which makes them sensitive to any kind of neutrino in¬ 
teraction and to all directions. The first hybrid analysis 
(“HI”) utilizes the outermost layer of optical modules as 
a veto against incoming muons, similar to the two shower 
searches SI and S2 described above. To ensure a suffi¬ 
ciently high veto probability, it is additionally required 
that more than 6,000 photo-electrons, integrated over all 
optical modules, be recorded. This value corresponds 
to a threshold of ~ 30 TeV in deposited energy. The 
second hybrid analysis (“H2”) employs a more complex 
veto algorithm that adapts to the energy of the event and 
additionally searches for early hits that could be caused 
by an incoming muon. As a result, the threshold for 
this event selection is significantly lower, ^ I TeV in 
deposited energy. A residual muon background contam¬ 
ination is present in the event samples of both searches, 
HI and H2. 
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Analysis HI uses data taken in 2010 — 2013 while anal¬ 
ysis H2 uses data taken in 2010 — 2012, both with the 
79-string and 86-string configurations of IceCube. For 
both event samples, energy and zenith angle information 
are available. Analysis H2 additionally provides data on 
the topology of the selected events; those in which more 
than 10 photo-electrons can be associated with an out¬ 
going muon track are classified as tracks, all others as 
showers.^ The atmospheric self-veto effect was imple¬ 
mented in the same way as for the shower searches (see 
above), based on Gaisser et al. (2014). 

2 . 4 . Dataset Overlap 

Data taken in one period is often subjected to multiple 
analyses, some of which may have overlapping event sam¬ 
ples. Specifically, among the analyses considered here, 
three (T2, HI, H2) use data taken between 2010 and 
2012 and have event samples that overlap. For the com¬ 
bined analysis presented here however, it is important 
that the different event samples be statistically indepen¬ 
dent, i.e. that individual events are not selected more 
than once. To achieve this, we identify the overlap be¬ 
tween the samples and remove the corresponding events 
from all but one of the samples, both in simulation and 
experimental data. In particular, we remove events that 
start inside the instrumented volume of IceCube from 
the event sample of analysis T2 and we consider only 
events from the event sample of analysis H2 for which 
fewer than 6,000 photo-electrons were recorded. With 
these adjustments, all event samples combined in the 
maximum-likelihood analysis are statistically indepen¬ 
dent. 

3. ANALYSIS METHOD 

In order to obtain a global picture of the astrophys- 
ical neutrino flux measured with IceCube, we analyze 
data from different event selections in a single maximum- 
likelihood analysis. From each of the searches for as- 
trophysical neutrinos introduced in the previous section, 
we use the final selection level sample of experimental 
data events as well as corresponding samples of simu¬ 
lated neutrino events and, if applicable, muon events. 
The data are binned in the available observables and an¬ 
alyzed together. For each component contributing to the 
measured flux, probability density functions (PDFs) are 
obtained from the simulated event samples by weight¬ 
ing them to different models. The sum of all PDFs is 
then compared to the distribution of experimental data. 
In the maximum-likelihood procedure, the parameters 
of the models are varied until the best agreement with 
the data is achieved. The models for the different back¬ 
ground components and the astrophysical component are 
described in sections 3.1 and 3.2, followed by a descrip¬ 
tion of the likelihood method (section 3.3). 

3 . 1 . Modeling of Background Components 

The relevant background components to all searches 
considered here are atmospheric neutrinos and atmo¬ 
spheric muons. We use the same models to describe the 
neutrino background in all of the different searches. The 
muon background is more specific and is hence modeled 
individually for each analysis. 

® This observable identifies >^60% (-^80%) of all charged-current 
i.'^-interactions as tracks at 10 TeV (1 PeV). 


Table 2 

Neutrino Background Models. 


Component 

Model 

Parameters 

Atmospheric conventional u 

1,3 

(pCOTiV 

Atmospheric prompt v 

2,3 

^prompt 


References. — (1) Honda et al. (2007); 
(2) Enberg et al. (2008); (3) Gaisser (2012) 


3 . 1 . 1 . Atmospheric Neutrinos 

The models and parameters for the atmospheric neu¬ 
trino components are summarized in Table 2. For the 
conventional atmospheric neutrinos, we use the HKKMS 
model by Honda et al. (2007) and for the prompt at¬ 
mospheric neutrinos, we use the ERS model by En¬ 
berg et al. (2008). Both models were slightly modi¬ 
fied (see Aartsen et al. 2014c) in order to be in accor¬ 
dance with measurements of the cosmic-ray spectrum in 
the energy range above 1 TeV; these modifications are 
based on the cosmic-ray composition model introduced 
in Gaisser (2012). The normalizations of both neutrino 
spectra (^iconv, ^prompt) are free parameters in the fit 
procedure.® 

3 . 1 . 2 . Atmospheric Muons 

The contamination of atmospheric muon background 
in the event samples of the track-based analyses T1 and 
T2 is negligible (Aartsen et al. 2014c, 2015a). In order 
to estimate the residual muon background in the event 
samples of the other analyses, simulations of cosmic ray 
air showers with CORSIKA (Heck et al. 1998) were car¬ 
ried out. For analysis SI, the expected muon contamina¬ 
tion amounts to ~57% of the observed number of events 
in the low-energy sample (Sla) and ~2% in the high- 
energy sample (Sib) (Aartsen et al. 2014d); for analysis 
S2 the expected contamination is ^12% (Schonwald et al. 
2013). The contributions at high energies were deter¬ 
mined by an extrapolation of the distribution at lower 
energies because of insufficient statistics at high energy. 
The magnitude of the muon background in the event 
sample of analysis HI was determined from experimen¬ 
tal data and is expected to be ~23% of the observed 
number of events (Aartsen et al. 2013, 2014e). Finally, 
the expected contamination in the analysis H2 is ~4% 
(Aartsen et al. 2015b). In all analyses, the estimation of 
the muon background is associated with large uncertain¬ 
ties. Accordingly, the normalizations of the muon back¬ 
ground contributions in each analysis are implemented 
as nuisance parameters, see section 3.3.2. 

3 . 2 . Modeling of the Astrophysical Component 

So far, the sources of the high-energy neutrino flux 
measured by IceCube have escaped identification (Aart¬ 
sen et al. 2014e,f). Here, we therefore test different mod¬ 
els for the astrophysical component, based on general 
assumptions about the neutrino flux. These models are 
listed in Table 3. 


® An update of the ERS model was recently presented by Bhat- 
tacharya et al. (2015), who find a lower prompt flux estimate. The 
spectral shape of the modified ERS model used here agrees well 
with that of the new model though, so that, because the absolute 
normalization is a free fit parameter, it can still be used without 
being in contradiction with the new model. 
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Table 3 

Models for the Astrophysical Neutrino Flux. 


Model 

Parameters 

single power law 

</>, 7 

differential 

4>1 - <t>9 

north-south 

<t>S, 7Si 4’n, In 

2-flavor 

06) 0/i-^T) 7 

3-flavor 

06) 0//) 0r) 7 


The simplest model is the “single power law” model. 
We assume that the astrophysical flux arrives isotropi¬ 
cally from all directions, that its flavor-ratio at Earth is 
= 1 : 1 : 1 , and that it can be described by a 
simple power law of the form 

(lOOTev) ■ 

<l>y denotes the all-flavor neutrino flux, (j) its value at 
100 TeV and 7 the power law spectral index. 

The “differential” model is based on the same assump¬ 
tions about isotropy and flavor composition, but models 
the astrophysical flux with nine independent basis func¬ 
tions, defined in nine logarithmically spaced energy in¬ 
tervals between 10 TeV and 10 PeV. The normalizations 
(^1 — (j)g of the basis functions are free fit parameters, 
while the energy spectrum in each interval is assnmed 
to be cx E~‘^. This model is similar to the procedure in 
Aartsen et al. (2014e, 2015b). 

In the “north-south” model, we relax the assumption 
of isotropy and allow for two independent astrophysical 
neutrino fluxes, one from the northern and one from the 
southern hemisphere (separated at declination 8 — 0 °), 
both following a spectrum as defined in equation ( 1 ). 
While this scenario might lack a good astrophysical mo¬ 
tivation, it does allow us to describe the flux separately in 
two hemispheres that are affected by different detector- 
related systematic effects. We refrain from testing more 
complex anisotropic models here because, having been 
selected for diffuse searches, the event samples used in 
this analysis currently do not contain unblinded right- 
ascension information. We do however note that our 
simple north-south model could be sensitive to certain 
anisotropic scenarios, like e.g. an additional astrophysi¬ 
cal component from the inner Galaxy. 

Finally, the assumptions about the flavor composition 
are weakened in the “2-flavor” and “3-flavor” models. In 
the 2 -flavor model, the and Vr dux are assumed to 
be equal at Earth, while the Ve dux is allowed to de¬ 
viate. This relation is, in good approximation, true for 
any flavor composition at the source if standard neutrino 
oscillations transform the neutrino dux during propaga¬ 
tion.' In the 3-davor model, the normalizations of the 
duxes of all three neutrino flavors are free parameters, al¬ 
lowing us to test for non-standard oscillation scenarios. 
In both the 2-flavor and 3-davor models, the duxes of 
the individual neutrino flavors are assumed to have the 
same energy dependence, i.e. oc E~^, where 7 is a free 
parameter. 


^ It is exactly true in the case of tribimaximal mixing, and valid 
up to ~20% for more realistic oscillation parameters. 


3 . 3 . The Maximum-likelihood Method 

In the maximum-likelihood method, the agreement be¬ 
tween experimental data and the simulated PDFs is es¬ 
tablished by means of the test statistic described in sec¬ 
tion 3.3.1. Section 3.3.2 lists the systematic uncertainties 
that we account for in the method. Finally, the calcu¬ 
lation of p-values for likelihood ratio and goodness-of-dt 
tests is explained in sections 3.3.3 and 3.3.4, respectively. 

3 . 3 . 1 . Test Statistic 

We use a binned Poisson-likelihood test statistic to 
compare the experimental data with the model predic¬ 
tions. The general dednition for the test statistic is 

-21nL = -2^1n[^^) , (2) 

where rii and Vi denote the number of observed and pre¬ 
dicted events in bin i, respectively. 

Systematic effects might change the PDFs that are dt 
to the data and could thus distort the results of the 
Ht. To avoid this, we have parametrized the impact on 
the PDFs of all relevant systematic effects and included 
them in the dt procedure as nuisance parameters. For 
each nuisance parameter 9 there is a prior, an additional, 
Gaussian-shaped penalty term in the likelihood function 
that penalizes deviations from the default central value 
0*. The width a[0] of the prior is based on the nncer- 
tainty associated with the systematic effect. With m 
nuisance parameters included, equation ( 2 ) now reads 


— 2 In L = —2 In 


m 


+ E 


\ ^[^3] 


■ ( 3 ) 


The individual systematic effects considered as nuisance 
parameters are discussed in greater detail in the next 
section. 


3 . 3 . 2 . Systematic Uncertainties 

The following systematic edects were included in the 
maximum-likelihood procedure (for a summary, see Ta¬ 
ble 4): 

— Cosmic-ray spectral index. Atmospheric neutri¬ 
nos are produced by cosmic rays hitting the atmo¬ 
sphere, hence their energy spectrum depends on the 
cosmic-ray energy spectrum. The uncertainty on 
the spectral index of this spectrum is implemented 
as a nuisance parameter that tilts the spectrum of 
atmospheric neutrinos by Ajci relative to the de¬ 
fault model. Note that positive values of Ay^r cor¬ 
respond to softer spectra. We use an uncertainty 
on the cosmic-ray spectral index of 0.05 (see e.g. 
Gaisser 2012). 

— Muon background normalization. The residual 
muon background was determined from simulations 
of cosmic-ray air showers and/or from data (see sec¬ 
tion 3.1.2). Due to the different methods used in 
the different individual analyses, systematic shifts 
can be uncorrelated between different analyses. We 
have therefore implemented one nuisance parame¬ 
ter for each search that varies the normalization 
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Table 4 

Systematic Uncertainties. 


Uncertainty 

Treatment 

Parameter 

Prior 

CR spectral index 

global 

Aqcr 

0 ± 0.05 

Atm. // background 

individual 


1 ± 0.5 

Energy scale 

individual 

Se 

1 ± 0.15 


Note. — The systematic uncertainties taken into 
account in this analysis. “Treatment” indicates whether 
this effect affects all analyses the same (“global”) or if 
the effect is treated different in each affected analysis 
(“individual”). 

of the residual muon background in that search 
only. Since the muon background is negligible for 
the track-based analyses T1 and T2, this amounts 
to a total of four muon background normalization 
parameters '/'m.hz)- The width 

of the prior is 50% for all of these parameters. 

— Energy scale shift. A shift of the energy scale could 
be introduced by several systematic effects. Be¬ 
cause the energy reconstruction is based on the 
number of detected photons, an uncertainty on the 
optical efficiency of the detector modules directly 
translates to an uncertainty on the energy scale. 
An imperfect ice model, which describes the scat¬ 
tering and absorption of photons in the ice and is 
used both in simulation and event reconstruction, 
could also lead to a shift in the energy scale. Since 
some analyses used older ice models than others, 
the effect could be different for the different analy¬ 
ses used here. Therefore, one energy scale nuisance 
parameter (j)E was implemented for each analysis, 
shifting the energy scale for that analysis only. We 
use a prior of 15% on each each energy scale pa¬ 
rameter. 

We also checked for the impact of a deviation of the 
conventional atmospheric electron-to-muon neutrino ra¬ 
tio from the model prediction, which would result from a 
mis-modeling of the kaon-to-pion ratio in atmospheric air 
showers. However, we found the impact on the parame¬ 
ters of the astrophysical neutrino flux to be negligible. 

3.3.3. Likelihood Ratio Tests 

Two different models HO and HI can be compared with 
likelihood ratio tests. For this, we use the quantity 

-2AlnL = -21nTHo/THi, (4) 

which Wilk’s theorem predicts to be x^-distributed with 
k degrees of freedom, where k is the difference in the 
number of parameters between HI and HO (Wilks 1938; 
Olive et al. 2014, chapter 38).* The distribution can de¬ 
viate from a x^"distribution if the sample size is small or 
a parameter is close to a physical bound. In this case, 
the exact distribution of —2 A In L can be computed from 
toy Monte-Carlo experiments, generated from the best- 
fit parameter values of model HO. The likelihood ratio 
test p-value is then given by the percentage of toy experi¬ 
ments that have a larger value of —2 A In L than observed 


® This requires the models to be nested, i.e. each point in the 
parameter space of HO can be accessed with the parameters of HI. 


in data. The likelihood ratio test p-values quoted in this 
article were obtained with this procedure. 

We also employ likelihood ratio tests to determine con¬ 
fidence intervals on the model parameters. For each pa¬ 
rameter value, we perform a likelihood ratio test between 
the model with the parameter constrained to this value 
and the model with the parameter unconstrained. This 
procedure is known as a profile likelihood scan. In the 
X^-approximation that we use here, the 68 % and 90% 
confidence level intervals are given by the values at which 
—2AlnL = 1 and — 2AlnL = 2.71, respectively. Simi¬ 
larly, 2 -dimensional confidence contours can be obtained 
by constraining two parameters simultaneously. 

3.3.4. Goodness-of-fit Tests 

With a slight modification, the quality of the fit of 
a single model can be assessed using the test statistic. 
Similar to equation (4), we define 

-21nL/L,at = -2^1n(e"--"-(i/,/n0”‘) , (5) 

i 

where Lgat is the likelihood of the model that exactly pre¬ 
dicts the observed outcome (i.e. Vi = Ui). After adding 
the same term as in equation (3), this quantity is mini¬ 
mized by the same parameter values that minimize the 
test statistic defined there. Moreover, in the large sample 
limit, the minimum value follows a x^"distribution and 
as such allows the calculation of a goodness-of-fit p-value 
(Olive et al. 2014, chapter 38). Since our samples are 
not large, we determine the distribution of the modified 
test statistic by generating toy Monte-Carlo experiments 
from the best-fit parameter values. Comparing this dis¬ 
tribution with the value observed in experimental data 
gives the goodness-of-fit p-value. 

4. RESULTS 

In this section, we present the results of the maximum- 
likelihood analysis for the different models of the astro- 
physical neutrino flux given in section 3.2 and Table 3. 

4.1. Single Power Law Model 

The best-fit parameter values for the single power law 
model are listed in Table 5, including all nuisance param¬ 
eters. Specifically, the all-flavor astrophysical neutrino 
flux is 

(j) = (6.7l(;^) • IQ-i® GeV^s-^sr-^cm-2 (6) 

at 100 TeV and the best-fit power law has a spectral 
index of 

7 = 2.50 ±0.09 . (7) 

This measurement is valid for neutrino energies between 
25 TeV to 2.8 PeV. This energy range was determined 
by successively removing events, ordered in energy, from 
the simulation data and repeating the fit with cj) and 
7 constrained to their best-fit values; its bounds de¬ 
note the energies at which the test statistic worsens by 
-2AlnL = 1. 

We obtain a reasonable p-value of 37.6% from a 
goodness-of-fit test for this model. A power law with 
a fixed index of 7 = 2 is disfavored with a significance 
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Table 5 

Best-Fit Parameter Values for the Single Power Law Model. 


Parana. 

Best fit 

GO 

3.L. 

90% 

C.L. 

Pull 

0conv 

1.10 

0.94 - 

1.31 

0.87- 

- 1.49 

_ 

^prompt 

0.00 

0.00 - 

1.04 

0.00- 

- 2.11 

— 

S 

6.7 

5.5 - 

7.8 

4.6 - 

- 8.6 

— 

7 

2.50 

2.41 - 

2.59 

2.35 - 

- 2.65 

- 

^7cr 

0.017 

-0.008 - 

- 0.041 

-0.023 

- 0.057 

0.34 


1.09 

0.72 - 

1.51 

0.52 - 

- 1.80 

0.18 


0.84 

0.31 - 

1.37 

0.00- 

- 1.71 

-0.32 

0H,H1 

1.12 

0.75 - 

1.54 

0.56 - 

- 1.84 

0.23 

0/j,H2 

1.27 

0.94 - 

1.61 

0.73 - 

- 1.84 

0.54 

<t>E,SI 

0.95 

0.88 - 

1.04 

0.84- 

- 1.12 

-0.34 

4>E,S2 

1.00 

0.88 - 

1.22 

0.83 - 

- 1.32 

0.03 

<i>E,T1 

1.02 

0.95 - 

1.09 

0.90- 

- 1.14 

0.10 

4>E,T2 

1.05 

0.97 - 

1.12 

0.93 - 

- 1.17 

0.30 

(pE,!!! 

0.96 

0.88 - 

1.06 

0.84- 

- 1.12 

-0.29 

4>E,112 

0.95 

0.86 - 

1.04 

0.81 - 

- 1.10 

-0.35 

Note. 

— ^ is the 

value 

of the 

all-flavor 

neutrino 

flux at 


100 TeV and is given in units of 10~^^ GeV“^s“^sr~^cm~^. 
0 conv and ^prompt are given as multiples of the model predic¬ 
tions (see Table 2). “Pull” denotes the deviation of a nuisance 
parameter from its default value in units of the prior width g. 


of 3.8 a {p = 0.0066%) in a likelihood ratio test with 
respect to the model with a free spectral index. We 
also tested a single power law model with a high-energy 
exponential cut-off as well as a model that consists of 
two isotropic astrophysical components, each described 
by a power law. Neither model gave a better descrip¬ 
tion of our data. A power law with a fixed index of 
7 = 2 and a high-energy exponential cut-off is still dis¬ 
favored with a signihcance of 2.1 a {p = 1.7%) with 
respect to the model with a free spectral index. The 
best-fit cut-off energy is (l.bl'lj!);®) PeV for a fixed spec¬ 
tral index of 7 = 2 , with a corresponding normalization 
of 4> = ( 5 . 2 ^i[' 5 ) • 10“^® GeV“^s“^sr“^cm“^. No cut-off 
is fitted for the best-fit spectral index of 7 = 2.5. 

The likelihood analysis favors a prompt atmospheric 
component equal to zero, with a 90% C.L. upper limit 
of 2.1 times the (modified) model prediction by En- 
berg et al. (2008). This limit is slightly higher than that 
obtained in Aartsen et al. (2015b); this is due partly to 
the way the data samples are combined here (cf. sec¬ 
tion 2.4) and to the different treatment of the energy 
scale uncertainty. 

The correlation between the parameters </> and 7 is vi¬ 
sualized in Figure 1. The behavior of these parameters 
as a function of the magnitude of the prompt component 
is shown in Figure 2. 

Figure 3 shows the results for the parameters ((>prompt, 
(j), and 7 again (row “combined”), together with results 
obtained from the application of the maximum-likelihood 
method described in this paper to the individual event 
samples. For the individual fits, the event samples that 
were reduced to remove overlap (see section 2.4) were 
restored to their original size. Furthermore, the nor¬ 
malization of the conventional atmospheric neutrino flux 
was treated as a nuisance parameter with a prior of 25% 
around the model prediction for the fits on the event 
samples of analyses SI, S2, and HI, since this compo¬ 
nent is not well constrained by these samples alone. A 
large prompt component instead of an astrophysical com¬ 
ponent is preferred in the fit of analysis S2 because these 



Figure 1. Profile likelihood scans around the best fit of the single 
power law model. The large panel displays a two-dimensional 
scan of the normalization (p and the spectral index 7 of the 
astrophysical neutrino flux; one-dimensional scans are shown in 
the small panels. The best-fit point is marked with “x” in the 
large panel, the dashed line shows the conditional best-fit value of 
(j) for each value of 7 . 

two components are close to degenerate if the only ob¬ 
servable is the deposited energy and the astrophysical 
spectrum is steep. Note that some of the results obtained 
with the individual samples differ slightly, although not 
signihcantly, from the originally published results (see 
references in Table 1). For instance, a softer spectral in¬ 
dex than presented in Aartsen et al. (2015a) is obtained 
for the sample T2 because only the high-energy data is 
used here; this difference is well within the uncertainty 
on this parameter. The somewhat harder spectral in¬ 
dex and lower normalization for sample HI measured in 
Aartsen et al. (2014e) is a result of the reduced energy 
range (> 60 TeV) used in the spectral fit there. 

Experimental and simulated data, weighted to the 
best-fit result, are shown in Figure 4 for all event sam¬ 
ples included in this analysis. The data are projected 
onto observables that are used in the likelihood fit. For 
the hybrid analyses HI and H2, we show the energy dis¬ 
tribution in two different zenith angle bins. 

Finally, the best-ht spectra for atmospheric and as¬ 
trophysical neutrinos in the single power law model are 
shown in Figure 5. 

4 . 2 . Differential Model 

Table 6 lists the best-fit parameter values for the differ¬ 
ential model. (j)i through (fg are the normalizations of the 
corresponding basis functions, each assumed to follow a 
spectrum oc E~^ and defined in logarithmically spaced 
energy intervals between 10 TeV and 10 PeV. The result¬ 
ing differential spectrum is shown in Figure 6 , together 
with the single power law model (cf. previous section). 

4 . 3 . North-South Model 

The best-fit values for the parameters of the astrophys¬ 
ical neutrino flux, separated into northern and southern 
hemisphere, are listed in Table 7. Most notable, the best- 
fit spectral index in the northern sky is = S.O^q®, 
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Figure 2. Profile likelihood scans around the best fit of the sin¬ 
gle power law model. The two large panels show two-dimensional 
scans of the normalization of the prompt atmospheric neutrino flux 
'?5>prompt and the normalization <j) (upper panel) and spectral index 
7 (lower panel) of the astrophysical neutrino flux, respectively; one¬ 
dimensional scans are shown in the small panels. See Figure 1 for 
further description. 

while in the southern sky it is 7 s = 2.56 ± 0.12. This 
discrepancy with respect to the single power law model 
corresponds to a statistical p-value of 13% (1.1 a). 

4.4. 2-Flavor and 3-Flavor Model 

For the 2-flavor model, the best-fit values for the as¬ 
trophysical Ve and -I- Vr flux are listed in Table 8 . The 
constraints on the astrophysical spectral index, which is 
assumed to be the same for both flavor components, are 
identical to those obtained in the single power law model. 
Furthermore, the results on the flavor composition do not 
depend strongly on the value of this parameter. We mea¬ 
sure an electron neutrino fraction of the astrophysical 
neutrino flux of0.18±0.11 at Earth. Figure 7 compares 
this value to fractions expected for different composition 
scenarios at the sources of the astrophysical flux. 

Finally, Table 9 gives the best-fit values for the astro- 
physical fluxes in the 3-flavor model. As in the 2 -flavor 
model, we find that the results do not change significantly 
when varying the spectral index of the astrophysical flux 
within its uncertainties. Note that none of the analy¬ 
ses considered here positively identifies charged-current 
tau neutrino interactions, which, depending on the tau 
decay mode, have a signature very similar to that of 
charged-current electron neutrino (~83%) or muon neu- 


Table 6 

Best-Fit Parameter Values for the Differential 
Model. 


Parameter 

Best fit 

68% C.L. 

90% C.L. 


9.3 

1.7- 17.3 

0.0 - 22.7 

<t>2 

22.6 

17.0 - 28.5 

13.5 - 32.5 

<t>3 

5.6 

2.4 - 9.2 

0.5 - 11.6 

4>i 

3.2 

0.8 - 5.9 

0.0 - 7.9 

<t>5 

4.3 

2.0 - 7.0 

0.8 - 9.0 

<f>6 

0.0 

0.0 - 1.5 

0.0 - 3.5 

4>7 

6.9 

4.5 - 9.7 

3.1 - 11.9 

<i>S, 

0.0 

0.0 - 1.5 

0.0 - 3.8 

4>9 

0.0 

0.0 - 0.6 

0.0 - 1.5 


Note. — 01 — are the all-flavor normaliza¬ 
tions (in of the individual basis functions, 

defined in nine logarithmically spaced energy inter¬ 
vals between 10 TeV and 10 PeV. They are given 
in units of 10~® GeVs“^sr~^cm~^. 


Table 7 

Best-Fit Parameter Values for the North-South 
Model. 


Parameter 

Best fit 

68% C.L. 

90% C.L. 

0iV 

2.1 

0.5 - 5.0 

0.1 - 7.3 

In 

2.0 

1.6 - 2.3 

1.2 - 2.5 

4>S 

6.8 

5.3 - 8.4 

4.4 - 9.5 

IS 

2.56 

2.44 - 2.67 

2.36 - 2.75 

Note. — 

- 0JV and 0 c are the all-flavor neutrino 

fluxes at 100 TeV in the northern 

and southern 

sky, respectively; 7 iv and 75 are the corresponding 

spectral indices. The fluxes are given in units of 

10-1® GeV 

-1 —1 —1 —9 

s sr cm . 



Table 8 

Best-Fit Parameter Values for the 2-Flavor 
Model. 

Parameter 

Best fit 68% C.L. 

90% C.L. 

0e 

1.3 0.5-2.1 

0.0- 2.6 

0/J,-l-T 

5.6 4.4-6.9 

3.7- 7.8 


Note. — 0e and 0 ^+t are the Ue. and 
+ t'r flux at 100 TeV, respectively. Both are 
given in units of 10“^® GeV~^s~^sr~^cm“^. 


trino (^17%) interactions. For this reason, the Vr flux 
is easily conflated with the and components in the 
3-flavor model. Figure 8 shows a profile likelihood scan 
of the flavor composition as measured on Earth, again 
compared to ratios expected for different source composi¬ 
tion scenarios. Performing likelihood ratio tests between 
these points and the best-fit hypothesis, we find p-values 
of 55% (0:1: 0), 27% (1:2: 0), and 0.014% (1:0: 0), 
respectively. Hence, our data are compatible with a pure 
muon neutrino composition and the generic pion-decay 
composition at the source, but reject a composition con¬ 
sisting purely of electron neutrinos at the source with a 
significance of 3.6 cr. 

5. DISCUSSION 

We have presented results of a maximum-likelihood 
analysis that is based on the combination of event sam¬ 
ples that were selected by six different studies designed 
to measure an astrophysical neutrino flux with IceCube. 
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Figure 3. Results for the parameters i/)prompt (left), (f> (center), and 7 (right) in the single power law model. The last row (“combined”) 
shows the results of the combined analysis, while the other rows show the results obtained by applying the fit method from this work to 
the individual event samples. The best fit is marked by the black data point, the shaded area corresponds to the 68% C.L. interval. 


Table 9 

Best-Fit Parameter Values for the 3-Flavor 
Model. 


Parameter 

Best fit 

68% C.L. 

90% C.L. 


2.9 

1.4- 

3.6 

0.0-4.2 

</>/4 

3.0 

2.4- 

3.7 

2.1 - 4.2 

(pj- 

0.0 

0.0- 

2.3 

0.0- 5.0 


Note. — 06, and (j>T are the Us, and 
Vt flux at 100 TeV, respectively. All are given 
in units of 10“^® GeV“'^s“^sr“^cm“^. 


Assuming that the astrophysical flux is isotropic and has 
a flavor ratio of Ug : : u,- = 1 : 1 : 1 at Earth, we found 

that the combined data are well described by a power law 
with spectral index —2.50 ± 0.09 between neutrino ener¬ 
gies of 25 TeV and 2.8 PeV. These are the most precise 
constraints on the spectrum of astrophysical neutrinos 
obtained so far. With present statistics, we found no 
evidence that more complex spectral shapes or features, 
such as a high-energy cut-off or multiple spectral compo¬ 
nents, are required to describe the data. However, astro- 
physically relevant deviations from a single power law are 
also not ruled out by the current data. Under the same 
assumptions of isotropy and flavor equality, we exclude a 
spectrum oc E~‘^, which is the benchmark prediction of 
the Fermi shock acceleration model (e.g. Gaisser 1990), 
with a significance of 3.8 cr {p = 0.0066%). Correspond¬ 
ingly, a spectrum oc E~^ with an additional high-energy 
exponential cut-off is disfavored with a significance of 
2.1 CT (p= 1.7%). 

Murase et al. (2013) pointed out that a soft spectrum 
(7 > 2.1 — 2 . 2 ) of the astrophysical neutrino flux might be 
difficult to accommodate with the observed extragalac- 
tic gamma-ray background (Ackermann et al. 2015), if 
the entire neutrino flux is produced in pp-interactions 
in extragalactic sources that are transparent to gamma 
rays and distributed according to the star-formation rate. 
The soft spectrum obtained in this analysis could there¬ 
fore be a first indication that at least part of the flux 
originates from other source classes, or that its spectral 
shape is more complex than assumed here. 

While the analysis favors the absence of a prompt at¬ 
mospheric neutrino component, a contribution to the flux 
at the level of the prediction by Enberg et al. (2008) is 
still allowed and would not alter the spectrum of the as¬ 
trophysical component significantly, as demonstrated in 
Figure 2. Figure 3 shows that the results are compatible 
with those found in the individual studies; differences to 
the originally published results are mainly due to differ¬ 
ent energy ranges used in the analysis. 


The strength of the astrophysical signal in different 
energy intervals is shown in the differential spectrum in 
Figure 6 . This spectrum suggests that it is mostly events 
with energies around 30 TeV that are responsible for the 
soft spectrum obtained in the analysis here. In fact, a 
previous analysis (Aartsen et al. 2014e) that did not in¬ 
clude data at these energies yielded a harder spectral 
index of —2.3 ± 0.3, but with larger uncertainties. The 
result is compatible with the one obtained here.® 

We have tested the hypothesis of isotropy by fitting a 
model with two astrophysical components, one in the 
northern and one in the southern sky. Compared to 
the all-sky result, the fit prefers a harder spectrum 
^ ■'^- 0 . 4 . 1 ) in the northern sky and a slightly softer 
spectrum (i<;-2.56±o.i2^ southern sky with a signif¬ 

icance of 1.1 a {p = 13%). The result is not conclusive; 
the discrepancy could be caused by a statistical fluctua¬ 
tion or by an additional component that is present in only 
one of the hemispheres (either an unmodeled background 
component, or e.g. a component from the inner Galaxy, 
although a single point source of the required strength 
to create the anisotropy anywhere in that region has al¬ 
ready been excluded (Adrian-Martinez et al. 2014)). Fur¬ 
ther analysis including right-ascension information will 
be helpful in testing the hypothesis of isotropy in the 
future. 

Finally, we performed a measurement of the flavor 
composition of the astrophysical neutrino flux. In a 
first test, we have measured the electron neutrino frac¬ 
tion at Earth in a tribimaximal mixing scenario, with 
equal and Vt fluxes at Earth. The best-fit fraction 
is 0.18 ± 0 . 11 , a value compatible with the fractions ex¬ 
pected from pion-decay sources (0.33) and muon-damped 
sources ( 0 . 22 ), but incompatible with that expected from 
neutron-beam sources (0.56), see Figure 7. In a sec¬ 
ond, more general test, we allow the normalizations of 
all three flavor components to vary independently and 
compare the result to compositions expected for differ¬ 
ent astrophysical scenarios in Figure 8 . In agreement 
with the first test, we find that pion-decay sources and 
muon-damped sources are well compatible with our data, 
while neutron-beam sources are disfavored with a signif¬ 
icance of 3.6 a {p = 0.014%). We do not find indications 
for non-standard oscillation scenarios. 

Previous measurements of the flavor composition 
were presented by Mena et al. (2014) and Palomares- 

® We have established the compatibility in a separate fit with¬ 
out the corresponding data set, i.e. without sample HI. The 
68% uncertainty interval for the spectral index obtained in this 
fit (—2.45 ± 0.10) overlaps with that obtained in Aartsen et al. 
(2014e). 
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Figure 4. Experimental data and model predictions. The simulated data are weighted to the best-fit result of the single power law model. 
The event sample ID is indicated in the bottom right corner of each panel (cf. Table 1). For analysis SI, both event samples Sla and Sib 
are combined in one panel (top center). For analysis H2, the event sample is shown separately for showers (H2a) and tracks {H2b). If 
additional cuts were applied, they are indicated in the top right corner of the panel. The lines show the conventional atmospheric neutrino 
flux (blue), the astrophysical neutrino flux (red), and the atmospheric muon background (yellow). The 90% C.L. upper limit on the prompt 
atmospheric neutrino flux is denoted by the dashed green line. The sum of all non-zero components is shown in gray. The black data points 
represent the experimental data, the errors bars denote 68% C.L. intervals as defined in Feldman ^ Cousins (1998). 
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E, [GeV] 


Conv. atmospheric {u^ + 1 ^ ) 

Prompt atmospheric {u^ +^'u, 90% C.L.) 
Astrophysica! {u^ + 1 / +v^) 


Figure 5. Best-fit neutrino spectra for the single power law model 
(all flavors combined). The blue and red shaded areas correspond 
to 68% C.L. allowed regions for the conventional atmospheric and 
astrophysical neutrino flux, respectively. The prompt atmospheric 
flux is fitted to zero, we show the 90% C.L. upper limit on this 
component instead (green line). 



Figure 6. Best-fit astrophysical neutrino spectra (all flavors com¬ 
bined). The red shaded area corresponds to the 68% C.L. allowed 
region for the single power law model (cf. Figure 5). The black 
data points show the result of the differential model; the horizontal 
bars denote the bin width, the vertical error bars denote 68% C.L. 
intervals. 


at source 
:0 1 : 0:0 
\ 


0 0.1 0.2 0.3 0.4 0.5 0.6 

z/g fraction at Earth 

Figure 7. Electron neutrino fraction measured at Earth in the 2- 
flavor model. The black point denotes the best-fit value, the filled 
bands show the 68% (green) and 90% (red) C.L. intervals. The 
dashed lines mark electron neutrino fractions expected for different 
flavor compositions at the source, assuming tribimaximal neutrino 
mixing angles. 
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Figure 8. Profile likelihood scan of the flavor composition 
at Earth. Each point in the triangle corresponds to a ratio 
Ve. ^II '■ as measured on Earth, the individual contribu¬ 

tions are read off the three sides of the triangle. The best-fit 
composition is marked with “x”, 68% and 95% confidence 
regions are indicated. The ratios corresponding to three flavor 
composition scenarios at the sources of the neutrinos, computed 
using the oscillation parameters in Gonzalez-Garcia et al. (2014, 
inverted hierarchy), are marked by the square (0:1:0), 
circle (1:2:0), and triangle (1 : 0 : 0), respectively. The 
best-fit composition obtained in an earlier IceCube analysis of 
the flavor composition (Aartsen et al. 2015c) is marked with a “-h”. 


Ruiz et al. (2015) (based on event sample HI, presented 
in Aartsen et al. 2014e), and by Palladino et al. (2015), 
Pagliaroli et al. (2015), and Aartsen et al. (2015c) (based 
on event samples that were extended with respect to HI, 
respectively). With respect to these measurements, the 
constraints presented here are significantly improved; we 
attribute this to the fact that the combined event sam¬ 
ple analyzed here contains a significant number of shower 
events as well as track events. Though the best-fit flavor 
composition obtained in Aartsen et al. (2015c) (white 
“-I-” in Figure 8) lies outside the 95% C.L. region, the 
68% C.L. region obtained here is completely contained 
within that obtained in the previous work, demonstrat¬ 
ing the compatibility of the two results. Because neither 
analysis was designed to identify tau neutrinos, a degen¬ 
eracy with respect to the ^T-fraction is observed in both, 
the slight preference towards a smaller ^^--contribution 
found here is likely connected to the slight differences in 
the energy distributions of the three neutrino flavors. In 
future, the identification of tau neutrinos will enable us 
to place stronger constraints on the flavor composition 
of the astrophysical neutrino flux. 
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APPENDIX A 

TABLE OF INTERACTION TYPES 

Table 10 lists the fraction of neutrino interaction types 
that contribute to the event samples introduced in sec¬ 
tion 2. 

APPENDIX B 

CHANGES MADE TO THE EVENT SELECTION OF 
SAMPLE S2 

The event sample of search S2 has been extended 
with respect to what was presented in Schonwald et al. 
(2013). Specifically, while only events with deposited 
energies larger than 38 TeV were analyzed in the origi¬ 
nal work, here we use all events with deposited energies 
above 7 TeV. To convince ourselves that this is possible, 
we have carefully evaluated the behavior of the back¬ 
ground components in the energy range that was added, 
in particular that of the residual atmospheric muon back¬ 
ground. Figure 9 shows the distribution of deposited 
energies at next-to-final selection level, where the event 
sample is still dominated by atmospheric muons. Above 
7 TeV, our simulation matches the observed data well. 
Below that energy, the selection efficiency rapidly drops 
and we observe a disagreement between simulation and 
experimental data. We conclude that the simulation data 
below 7 TeV should not be used for analysis. 

The distribution of deposited energies at final selection 
level is shown in Figure 10, the data above 7 TeV is used 
in this work. Note that unlike in Figure 4, the simulation 
data has not been fitted to the experimental data for this 
Figure, but shows baseline predictions. As before, pre¬ 
dictions for atmospheric neutrinos are from Honda et al. 
(2007) (conventional) and Enberg et al. (2008) (prompt), 
both updated to the cosmic-ray composition model by 
Gaisser (2012) and corrected for the self-veto effect us¬ 
ing the calculation from Gaisser et al. (2014). The resid¬ 
ual muon background has been determined with a COR- 
SIKA simulation (Heck et al. 1998), the simulated events 
remaining in the final sample are indicated by the pur¬ 
ple histogram. As this simulation data lacks statistics 
in particular at high energies, we have parametrized and 
extrapolated the distribution (yellow line). This extrap¬ 
olation is used as a pdf in the combined analysis pre¬ 
sented here. We note that at final level, the residual 



Figure 9. Distribution of deposited energy for the event sample of 
search S2 at next-to-final selection level. The lines show predictions 
for the conventional atmospheric neutrino flux (blue), the prompt 
atmospheric neutrino flux (green) and atmospheric muons (purple). 
The black data points represent experimental data, the error bars 
denote 68% C.L. intervals as defined in Feldman &: Cousins (1998). 



Figure 10. Same as Figure 9, but at final selection level. Ad¬ 
ditionally, an extrapolation of the muon background is shown in 
yellow. The gray filled histogram shows the sum of the two at¬ 
mospheric neutrino components and the extrapolated muon back¬ 
ground. The original work used events above 38 TeV (black arrow), 
here we use the data above 7 TeV (red arrow). 

muon background is never the dominant component in 
the entire energy range and conclude that it is possible 
to use the data in the range between 7 and 38 TeV. 

Figure 10 also shows that even though the analysis 
in Schonwald et al. (2013) did not find an astrophysi- 
cal component in the neutrino flux, there is an excess of 
experimental data above the baseline background predic¬ 
tions. This excess was attributed to a very large prompt 
component in the original work, although with very large 
uncertainties. We qualitatively reproduce this result in 
our individual fit of sample S2 (see Figure 3, fourth row); 
we attribute the quantitative differences to the different 
energy range analyzed and the fact that the self-veto ef¬ 
fect was neglected in the original work. Furthermore, 
we observe that the excess is also compatible with an 
astrophysical flux as measured in the combined analysis 
presented here (cf. Figure 4, panel S2). This shows that 
when analyzed individually, the data of sample S2 can¬ 
not be used to distinguish between a large prompt flux 
and a steep astrophysical flux. 
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Table 10 

Fraction of Interaction Types per Event Sample 


10 TeV <Eu < 100 TeV 100 TeV <E^ < I PeV 

Ue Vr Ve. 

CC NC CC NC CC NC GR CC NC CC NC CC NC GR 


Conventional Neutrinos 


T1 

0.0 

0.0 

100.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

100.0 

0.0 

0.0 

0.0 

0.0 

T2a 

0.0 

0.0 

100.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

100.0 

0.0 

0.0 

0.0 

0.0 

SI 

7.7 

1.5 

48.8 

42.0 

0.0 

0.0 

0.0 

5.2 

1.5 

40.8 

52.5 

0.0 

0.0 

0.0 

S2 

17.0 

1.6 

32.7 

48.7 

0.0 

0.0 

0.0 

8.8 

2.3 

30.2 

58.7 

0.0 

0.0 

0.0 

HI 

15.4 

0.6 

65.7 

18.3 

0.0 

0.0 

0.0 

6.9 

0.7 

73.1 

19.3 

0.0 

0.0 

0.0 

H2a'^ 

13.6 

1.9 

30.5 

54.0 

0.0 

0.0 

0.0 

2.4 

2.5 

23.2 

71.9 

0.0 

0.0 

0.0 

H2b=- 

0.0 

0.0 

100.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

99.2 

0.8 

0.0 

0.0 

0.0 


Astrophysical Neutrinos {E 


T1 

0.0 

0.0 

96.5 

0.0 

3.5 

0.0 

0.0 

0.0 

0.0 

91.7 

0.0 

8.3 

0.0 

0.0 

T2® 

0.0 

0.0 

97.4 

0.0 

2.6 

0.0 

0.0 

0.0 

0.0 

96.9 

0.0 

3.1 

0.0 

0.0 

SI 

43.1 

9.0 

8.0 

7.0 

24.7 

8.1 

0.1 

34.2 

10.3 

7.8 

11.0 

26.7 

9.7 

0.3 

S2 

46.3 

5.9 

3.9 

5.9 

31.9 

6.0 

0.1 

35.9 

9.7 

4.7 

9.4 

29.6 

10.4 

0.3 

HI 

58.1 

2.8 

9.9 

2.7 

23.6 

2.7 

0.2 

37.8 

4.9 

19.4 

5.0 

27.4 

5.1 

0.4 

H2a^ 

41.0 

6.6 

4.5 

6.5 

33.7 

7.5 

0.2 

12.2 

15.2 

5.0 

15.1 

34.0 

18.4 

0.1 

H2b=- 

0.1 

0.0 

90.7 

0.0 

9.1 

0.1 

0.0 

0.3 

0.1 

79.6 

0.8 

18.5 

0.7 

0.0 


Note. — Fraction, in %, of neutrino interaction types contributing to the individual event samples. CC = 
charged-current interaction; NC = neutral-current interaction; GR = Glashow resonance (Pe +e“ ^ W~). 
The top table shows numbers for a spectrum weighted to the conventional atmospheric neutrino flux of 
Honda et al. (2007); the bottom table is for an astrophysical power-law spectrum with index —2.5 and 
flavor composition i^e • = 1 • 1 • 1- 

^ Note that these event samples were adjusted to remove overlap between samples T2, HI, and H2, see 
section 2.4. 
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